#include "seedcounter.h"

SeedCounter::SeedCounter()
{}

SeedCounter::~SeedCounter()
{}

bool SeedCounter::addFile(const QString &file)
{
    total = 0;
    masksize = 3;
    red_thres = 250;
    green_thres = 250;
    blue_thres = 254;
    seed_thres = 191;
    min_area = 50;
    max_area = 3000;
    PI = 3.14159265358979323846;

    blobareaavg = 0;
    blobperimeteravg = 0;
    blobroundnessavg = 0.0;
    blobaspectavg = 0.0;
    blobmajaxisavg = 0;
    blobminaxisavg = 0;

    return isFormat(file);
}

/*
 *	To test if input file is of particular type in which both
 *	top and bottom part of petri dish are visible. Otherwise
 *	output error messages
 */
bool SeedCounter::isFormat(const QString &file)
{
    bool isFormat = false;

    if(file.length() <= 0) return false;

    Mat img = imread(file.toStdString());
    if(!img.data) return false;

    int d = 0;
    Mat img_resized = resizeImage(img);
    Mat img_mask = findROIMask(img_resized, d);

    if(d > img_mask.cols/2)
    {
        run(img_resized, img_mask);
        isFormat = true;
    }
    return isFormat;
}

void SeedCounter::run(const Mat& img_resized, const Mat& img_mask)
{
    Mat img_cropped;
    Mat img_filtered;
    Mat img_segmented;

    if(img_resized.data)
    {
        // crop the image using the image mask
        img_cropped = cropImage(img_resized, img_mask);
        // apply noise filter to the working image
        img_filtered = customNoiseFilter(img_cropped, red_thres, green_thres, blue_thres);
        // segment image
        img_segmented = automatedSegmentation(img_filtered, img_mask);
        resultImg = calculateSeeds(img_segmented, img_resized, min_area, max_area);
    }
    else
    {
        // Image was not loaded correctly
        qDebug(QString("Error loading image.").toStdString().c_str());
    }
}


/*
 *	Takes in an original color image.
 *	Returns a resized color image.
 */
Mat SeedCounter::resizeImage(const Mat &image)
{
    Mat _resized, result;

    resize(image, _resized, Size(), 0.5, 0.5, CV_INTER_AREA);
    result = Mat(_resized, Rect(0,0,(int)(_resized.cols/1.75), (int)(_resized.rows)));

    return result;
}
/*
 *	Takes in an original color image.
 *	Returns a binary mask of image with ROI as well diameter
 *	of located circle
 */
Mat SeedCounter::findROIMask(const Mat &image, int& d)
{
    Mat croppedcopy = image.clone();
    Mat croppedgray = Mat(image.size(), CV_8UC1);

    //Need to crop the original image to only capture the petri dish
    cvtColor(croppedcopy, croppedgray, CV_BGR2GRAY);
    GaussianBlur(croppedgray, croppedgray, Size(9, 9), 0.5, 0.5);
    vector<Vec3f> circles;
    HoughCircles(croppedgray, circles, CV_HOUGH_GRADIENT, 2, croppedgray.cols, 50, 100);

    int radius = 0;
    for( size_t i = 0; i < circles.size(); i++ )
    {
         Point center(cvRound(circles[i][0]), cvRound(circles[i][1]));
         radius = cvRound(circles[i][2]);
         // draw the circle outline
         circle(croppedcopy, center, radius, Scalar(255,255,255), -1);
    }

    // Diameter of located circle
    d = 2 * radius;
    cvtColor(croppedcopy, croppedgray, CV_BGR2GRAY);
    threshold(croppedgray, croppedgray, 254, 255, CV_THRESH_BINARY);

    return croppedgray;
}
/*
 *	Takes in an resized color image.
 *	Returns a cropped color image.
 */
Mat SeedCounter::cropImage(const Mat &image, const Mat &imageMask)
{
    Mat _croppedfinal = Mat::zeros(image.size(), CV_8UC3);
    image.copyTo(_croppedfinal, imageMask);
    return _croppedfinal;
}
/*
 *	Takes in a color image and three threshold values for red, green, and blue.
 *	Returns a noise filtered color image.
 */
Mat SeedCounter::customNoiseFilter(const Mat &image, int r_threshold, int g_threshold, int b_threshold)
{
    int iter = 10;
    int masksize = 3;

    vector<Mat> planes;

    Mat img_working = Mat(image.size(), CV_8UC3);
    Mat r_median = Mat(image.size(), CV_8UC1);
    Mat g_median = Mat(image.size(), CV_8UC1);
    Mat b_median = Mat(image.size(), CV_8UC1);
    Mat r_mask = Mat(image.size(), CV_8UC1);
    Mat g_mask = Mat(image.size(), CV_8UC1);
    Mat b_mask = Mat(image.size(), CV_8UC1);

    //need to split the color into three r, g, b color channels
    split(image, planes);

    //normalize the image using histogram equalization
    equalizeHist(planes[0], r_mask);
    equalizeHist(planes[1], g_mask);
    equalizeHist(planes[2], b_mask);

    //create binary image capturing only the noise
    threshold(r_mask, r_mask, r_threshold, 255, CV_THRESH_BINARY);  //working
    threshold(g_mask, g_mask, g_threshold, 255, CV_THRESH_BINARY);  //working
    threshold(b_mask, b_mask, b_threshold, 255, CV_THRESH_BINARY);  //working

    //apply median filter to be used as a source when applying copy image using the mask
    medianBlur(planes[0], r_median, masksize);
    medianBlur(planes[1], g_median, masksize);
    medianBlur(planes[2], b_median, masksize);

    //using iterations to fill in the noise with neighboring pixels
    for(int i = 0;i < iter;i++)
    {
        medianBlur(r_median, r_median, masksize);
        medianBlur(g_median, g_median, masksize);
        medianBlur(b_median, b_median, masksize);
    }

    //copy only the areas in the mask to the working image
    r_median.copyTo(planes[0], r_mask);
    g_median.copyTo(planes[1], g_mask);
    b_median.copyTo(planes[2], b_mask);

    //merge the filtered channels back into a color image
    merge(planes, img_working);

    //apply median filter to working image and apply segmentation
    medianBlur(img_working, img_working, masksize);

    return img_working;
}

/*
 *	Takes in a color image
 *	Returns a segmented binary image
 */
Mat SeedCounter::automatedSegmentation(const Mat &image, const Mat &imageMask)
{
    Mat img_gray = Mat(image.size(), CV_8UC1);
    Mat imageMaskCpy = imageMask.clone();

    cvtColor(image, img_gray, CV_BGR2GRAY);
    equalizeHist(img_gray, img_gray);
    adaptiveThreshold(img_gray, img_gray, 255, CV_ADAPTIVE_THRESH_MEAN_C, CV_THRESH_BINARY_INV, 39, 5);
    bitwise_not(imageMaskCpy, imageMaskCpy);
    imageMask.copyTo(img_gray, imageMaskCpy);

    Mat kern = getStructuringElement(MORPH_RECT, Size(7, 7), Point(3,3));
    morphologyEx(img_gray, img_gray, MORPH_OPEN, kern);

    return img_gray;
}

/*
 *	Takes in segmented binary image, color image, min and max area to be considered a blob of interest
 *	Returns a color image with seeds outlined in red (single seed) or blue (multiple seeds)
 */
Mat SeedCounter::calculateSeeds(const Mat &imageGray, const Mat &imageColor, int min_area, int max_area)
{
    CBlobResult blobs;

    Mat _imageSeeds = imageColor.clone();
    Mat _imageBlobs = imageGray.clone();

    IplImage imageSeeds = _imageSeeds;
    IplImage imageBlobs = _imageBlobs;

    blobs = CBlobResult(&imageBlobs, NULL, 0);
    blobs.Filter(blobs, B_EXCLUDE, CBlobGetArea(), B_LESS, min_area);
    blobs.Filter(blobs, B_EXCLUDE, CBlobGetArea(), B_GREATER, max_area);
    ////////////////printf("NUMBER of Blobs ARE: %d\n",blobs.GetNumBlobs());

    CBlob *blob;
    CvPoint p;

    int totalarea = 0;
    int currentarea = 0;
    int currentperimeter = 0;
    int avgarea = 0;
    int blobarea = 0;
    int blobperimeter = 0;
    int blobmajaxis = 0;
    int blobminaxis = 0;
    int boundheight = 0;
    int boundwidth = 0;
    int seedcount = 0;
    int color = 0;

    double blobroundness = 0.0;
    double blobaspect = 0.0;


    //First iteration to determine the total area of all the detected blobs
    for(int i=0;i<blobs.GetNumBlobs();i++) {
            blob = blobs.GetBlob(i);
            currentarea = (int)blob->Area();
            totalarea = totalarea+currentarea;
    }
    //Calculate the average area of the blobs
    if(blobs.GetNumBlobs()>0) {
            avgarea = totalarea/blobs.GetNumBlobs();
    }

    //Recalculate the seed count by counting only the blobs that has seed potential
    totalarea = 0;
    for(int i=0;i<blobs.GetNumBlobs();i++) {
            blob = blobs.GetBlob(i);
            currentarea = (int)blob->Area();
            if(currentarea >= avgarea*2) {
                    totalarea = totalarea+currentarea;
                    seedcount = seedcount + currentarea/avgarea;
            } else if (currentarea >= avgarea/2.5) {
                    totalarea = totalarea+currentarea;
                    seedcount++;
            }
    }

    total = seedcount;
    //Recalculate the average area of the seed
    if(blobs.GetNumBlobs()>0) {
            avgarea = totalarea/seedcount;
    }
    /////////////printf("NUMBER of Seeds ARE: %d\n", seedcount);

    //Second iteration to count the seeds applying the average area to the seed counting logic
    for(int i=0;i<blobs.GetNumBlobs();i++)
    {
        color = (i*5 > 255)? 255 : i*5;
        blob = blobs.GetBlob(i);

        CvBox2D box;
        box.angle = CBlobGetOrientation()(*blob);
        box.center = cvPoint2D32f(CBlobGetXCenter()(*blob), CBlobGetYCenter()(*blob));

        currentarea = (int)blob->Area();
        currentperimeter = (int)blob->Perimeter();
        boundheight = (int)CBlobGetMajorAxisLength()(*blob);
        boundwidth = (int)CBlobGetMinorAxisLength()(*blob);

        // Each blob's centroid point
        p = cvPoint(box.center.x, box.center.y);


        if(boundheight > boundwidth)
        {
            box.size.height = boundheight*2;
            box.size.width = boundwidth*2;
        }
        else
        {
            box.size.height = boundwidth*2;
            box.size.width = boundheight*2;
        }
        if(currentarea >= avgarea*2)
        {
            // Need to fill each seed with different color
            cvEllipseBox(&imageSeeds, box, CV_RGB(0, 0, 255), 10);
            //calculate blob data
            blobarea = currentarea;
            blobperimeter = currentperimeter;
            blobmajaxis = (int)box.size.height;
            blobminaxis = (int)box.size.width;
            blobroundness = pow((double)blobperimeter, 2.0) / (4.0 * PI * (double)blobarea);
            blobaspect = (double)blobmajaxis/(double)blobminaxis;
            //////////////////printf("Seed #%d: %d [area]; %d [perimeter]; %d [major]; %d [minor]; %.2f [roundness]; %.2f [aspect]\n", i, blobarea, blobperimeter, blobmajaxis, blobminaxis, blobroundness, blobaspect);
            blobareaavg = blobareaavg + blobarea;
            blobperimeteravg = blobperimeteravg + blobperimeter;
            blobmajaxisavg = blobmajaxisavg + blobmajaxis;
            blobminaxisavg = blobminaxisavg + blobminaxis;
            blobroundnessavg = blobroundnessavg + blobroundness;
            blobaspectavg = blobaspectavg + blobaspect;

        }
        else if (currentarea >= avgarea/2.5)
        {
            //map->insert(QString::number(i), *box);
            // need to fill each seed with different color
            cvEllipseBox(&imageSeeds, box, CV_RGB(255-color,255-i,color), 5);
            //calculate blob data
            blobarea = currentarea;
            blobperimeter = currentperimeter;
            blobmajaxis = (int)box.size.height;
            blobminaxis = (int)box.size.width;
            blobroundness = pow((double)blobperimeter, 2.0) / (4.0 * PI * (double)blobarea);
            blobaspect = (double)blobmajaxis/(double)blobminaxis;
            //////////////printf("Seed #%d: %d [area]; %d [perimeter]; %d [major]; %d [minor]; %.2f [roundness]; %.2f [aspect]\n", i, blobarea, blobperimeter, blobmajaxis, blobminaxis, blobroundness, blobaspect);
            blobareaavg = blobareaavg + blobarea;
            blobperimeteravg = blobperimeteravg + blobperimeter;
            blobmajaxisavg = blobmajaxisavg + blobmajaxis;
            blobminaxisavg = blobminaxisavg + blobminaxis;
            blobroundnessavg = blobroundnessavg + blobroundness;
            blobaspectavg = blobaspectavg + blobaspect;
        }
    }
    if(seedcount > 0)
    {
        //calculate average blob data
        blobareaavg = blobareaavg/seedcount;
        blobperimeteravg = blobperimeteravg/seedcount;
        blobmajaxisavg = blobmajaxisavg/seedcount;
        blobminaxisavg = blobminaxisavg/seedcount;
        blobroundnessavg = blobroundnessavg/seedcount;
        blobaspectavg = blobaspectavg/seedcount;
    }
    ///////////////printf("Averages: %d [area]; %d [perimeter]; %d [major]; %d [minor]; %.2f [roundness]; %.2f [aspect]\n", blobareaavg, blobperimeteravg, blobmajaxisavg, blobminaxisavg, blobroundnessavg, blobaspectavg);

    blobs.ClearBlobs();

    return _imageSeeds;
}

int SeedCounter::numberOfSeeds()
{
    return total;
}

int SeedCounter::avgArea()
{
    return blobareaavg;
}

int SeedCounter::avgPeri()
{
    return blobperimeteravg;
}

double SeedCounter::avgRoundness()
{
    return blobroundnessavg;
}

double SeedCounter::avgAspect()
{
    return blobaspectavg;
}

int SeedCounter::avgMajor()
{
    return blobmajaxisavg;
}

int SeedCounter::avgMinor()
{
    return blobminaxisavg;
}


